Tarea 5: MCMC

CC6104: Statistical Thinking

Integrantes :

Franco Cattani
Juan Molina

Cuerpo Docente:

Profesor: Felipe Bravo M.
Auxiliar: Martín Paredes, Benjamín Farías
Ayudantes: Scarleth Betancurt, Emilio Díaz, Kevin Iturra, Renzo Zanca

Índice:

Objetivo
Instrucciones
Referencias
Primera Parte: Preguntas Teóricas
Segunda Parte: Elaboración de Código

Objetivo

Bienvenid@s a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

La tarea se realiza en grupos de máximo 2 personas. Pero no existe problema si usted desea hacerla de forma individual.
La entrega es a través de u-cursos a más tardar el día estipulado en la misma plataforma. A las tareas atrasadas se les descontará un punto por día.
El formato de entrega es este mismo Rmarkdown y un html con la tarea desarrollada. Por favor compruebe que todas las celdas han sido ejecutadas en el archivo html.
Al momento de la revisión tu código será ejecutado. Por favor verifica que tu entrega no tenga errores de compilación.
No serán revisadas tareas desarrolladas en Python.
Está PROHIBIDO la copia o compartir las respuestas entre integrantes de diferentes grupos.
Pueden realizar consultas de la tarea a través de U-cursos y/o del canal de Discord del curso.

Referencias:

Slides de las clases:

Videos de las clases:

Introduction to Bayesian Inference: video1 video2 video3 video4 video5 video6 video7
Summarizing the Posterior: video1 video2 video3 video4

Documentación:

Pregunta 1: Model Evaluation and Information Criteria

Explique cómo cross-validation, criterios de información y regularización ayudan a mitigar o medir los problemas de underfitting y overfitting.

Respuesta

En el caso de cross-validation, lo hace mediante testear en pequeños subsets de la data e ir testeando con cada uno y entrenando con lo demás. Básicamente, es un proceso iterativo que testea en todos los datos de a poco, lo que permite al modelo no memorizar la data de entrenamiento ya que esta cambia en cada iteración.

Regularización evita el overfitting y underfitting mediante ajusta un valor lambda el cual se usa en la función likelihood, el cual penaliza para parámetros con valores grandes. Esto genera que los priors sean más escépticos según como se ajuste el parámetro (Si es menor, se evita underfitting, si es mayor, se evita overfitting).

Los criterios de información, en particular AIC, utilia el número de parámetros para penalizar a los modelos que tienen muchos parámetros (overfitting) o compensar para modelos que se ajustan mejor a los datos (underfitting). Esto lo hace mediante estimar su valor como la suma entre la likelihood y el doble del número de parámetros (lo cual permite el ajuste). Este valor equivale a una estimación de la out-of-sample deviance promedio.

Pregunta 2: Directed Graphical Models

Diseñe una DAG para un problema causal inventado por usted de al menos 5 nodos (puede basarse en el ejemplo de Eugene Charniak) usando Dagitty y considere que la DAG tenga al menos: una chain, un fork y un collider. Muestre con dagitty todas las independencias condicionales de su DAG. Explique las independencias usando las reglas de d-separación.

Respuesta

install.packages("dagitty", repos = "https://cloud.r-project.org/")

## package 'dagitty' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\juanm\AppData\Local\Temp\Rtmpg5DOfZ\downloaded_packages

Usaremos un ejemplo en un contexto médico. Una enfermedad puede producir fiebre o dolores de cabeza. Si el paciente está deshidratado, puede ser por fiebre o por náuseas. Esto junto con las independencias condicionales se visualiza en el siguiente código que utiliza dagitty:

library(dagitty)

## Warning: package 'dagitty' was built under R version 4.4.2

# s := sickness, f := fever, h := headache, d := dehydration, n := nausea
over <- dagitty("dag{ s -> f; s -> h; f -> d; n -> d }")

coordinates(over) <- list(x = c(s = 2, f = 3, h = 1, d = 2, n = 2),
                          y = c(s = 4, f = 3, h = 3, d = 2, n = 1))

plot(over)

cat("Las independencias condicionales son:\n")

## Las independencias condicionales son:

impliedConditionalIndependencies(over)

## d _||_ h | s
## d _||_ h | f
## d _||_ s | f
## f _||_ h | s
## f _||_ n
## h _||_ n
## n _||_ s

Luego explicamos cada independencias condicional.

d || h | s

Esta independencia nos indica que la deshidratación es independiente del dolor de cabeza, dado que conocemos que se tiene la enfermedad. El dolor de cabeza no influye en la deshidratación por lo que está independencia es evidente y por regla de d-separación debería estar d-separados condicionado por s, ya que s está en un fork en el camino entre d y h y es parte de la evidencia:

dseparated(over, "d", "h", c("s"))

## [1] TRUE

Lo mismo pasa con la siguiente, en este caso respecto a la fiebre f, y f aquí es parte de una cadena:

d || h | f

dseparated(over, "d", "h", c("f"))

## [1] TRUE

La siguiente nos muestra que la deshidratación es independiente de la enfermedad dado la fiebre. Si sabemos que tiene fiebre, entonces la enfermedad no impacta en la probabilidad de deshidratación. Esto significa que deberían estar d-separados.

d || s | f

dseparated(over, "d", "s", c("f"))

## [1] TRUE

La siguiente nos dice que la fiebre es independiente del dolor de cabeza dado que existe la enfermedad. Esto se explica ya que la fiebre y el dolor de cabeza son independientes si se sabe que se tiene la enfermedad. Si no, se podría saber que se tiene fiebre, lo cual podría significar que se tiene la enfermedad y esto abre a la posibilidad de que se tenga dolor de cabeza (Por ende no son independientes si no se sabe que se tiene la enfermedad). La enfermedad es un fork entre ellos dos, pero es parte de la evidencia, por lo que deberían estar d-separados:

f || h | s

dseparated(over, "f", "h", c("s"))

## [1] TRUE

Esta independencia nos dice que la fiebre es independiente de las nauseas. Esto tiene sentido ya que las nauseas solo causan deshidratación, pero esta última no esta involucrada en la relación entre la fiebre y las nauseas. Esto último nos indica que d no es parte de la evidencia y es un collider, por lo que deben estar d-separadas:

f || n

dseparated(over, "f", "n", c())

## [1] TRUE

La siguiente nos indica que el dolor de cabeza es independiente de las nauseas. No existe relación causal entre las nauses y el dolor de cabeza, por lo que esta independencia es clara. Entre estos nodos hay tanto un collider, un fork y una cadena. El collider d no esta en la evidencia, y los otros dos tampoco, por lo que deben estar d-separados:

h || n

dseparated(over, "h", "n", c())

## [1] TRUE

Por último, una independencia entre las nauseas y la enfermedad. Si se tiene la enfermedad, no nos dice nada sobre las nauseas. Entre n y ese esta el collider d, que no es parte de la evidencia y la cadena f que tampoco lo es, por lo que deben estar d-separados:

dseparated(over, "n", "s", c())

## [1] TRUE

n || s

Segunda Parte: Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
library(dplyr)
library(tidyr)

# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)

# Manipulación de varios plots en una imagen.
library(gridExtra)

# Análisis bayesiano
library("rethinking")

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')

Pregunta 1: MCMC

El objetivo de esta pregunta es lograr samplear, mediante la técnica de MCMC, la distribución gamma.

En general la distribución gamma se denota por \(\Gamma(\alpha,\beta)\) donde \(\alpha\) y \(\beta\) son parámetros positivos, a \(\alpha\) se le suele llamar “shape” y a \(\beta\) rate La densidad no normalizada de una distribución gamma esta dada por:

\[ f(x\mid \alpha,\beta) = \begin{cases} x^{\alpha -1}e^{-\beta x} ~ &\text{ si } x > 0\\ 0 ~&\text{si } x \leq 0 \end{cases} \] donde \(\Gamma(\alpha)\) es una constante, usualmente se le llama función gamma.

Implemente el algoritmo de metropolis hasting utilizando \(q(\theta^{(t)} \mid \theta^{(t-1)}) = \mathcal{N}(\theta^{t-1},1)\), importante su función tiene que poder recibir:
- La condición inicial \(\theta_{0}\).
- Cantidad de repeticiones.
- \(\alpha\) y \(\beta\).
Escriba el algoritmo sin utilizar implementenaciones de la distribución gamma en r.

De ahora en adelante considere \(\alpha = 5\) y \(\beta = \frac{1}{5}\).

Considere \(\theta_{0} = 1\), grafique el histograma que obtiene para distintas cantidad de repeticiones, entre sus pruebas tiene que tener al menos una que tenga del orden de \(10^5\) repeticiones ¿Cómo cambia la distribución en funcion de las repeticiones?
Considere distintos valores de \(\theta_{0}\) y una cantidad de repeticiones grande (del orden de \(10^5\)), grafique las distribuciones que obtenga comente sus resultados ¿Es importante la condición inicial del algoritmo?.
Utilizando \(\theta_{0}\) y cantidad de repeticiones conveniente (de acuerdo a lo que obtuvo en las partes anteriores) compare con la distribución real. Obs: En esta parte puede utilizar la distribución gamma que tiene implementado R para comparar con lo que usted realizó.

Respuesta

# Algoritmo Metropolis-Hastings
metropolis_hastings <- function(theta0, alpha, beta, iterations) {
  gamma_density <- function(x, alpha, beta) {
    if (x > 0) {
      return(x^(alpha - 1) * exp(-beta * x))
    } else {
      return(0)
    }
  }
  
  samples <- numeric(iterations)
  samples[1] <- theta0
  
  for (t in 2:iterations) {
    theta_star <- rnorm(1, mean = samples[t - 1], sd = 1)
        if (theta_star > 0) { 
      acceptance_ratio <- gamma_density(theta_star, alpha, beta) / 
                          gamma_density(samples[t - 1], alpha, beta)
    } else {
      acceptance_ratio <- 0
    }
    
    u <- runif(1)
    if (u < acceptance_ratio) {
      samples[t] <- theta_star
    } else {
      samples[t] <- samples[t - 1]
    }
  }
  
  return(samples)
}

# Configuración inicial
alpha <- 5
beta <- 15 
theta0 <- 1

# Diferentes cantidades de repeticiones
samples_1k <- metropolis_hastings(theta0, alpha, beta, 1000)
samples_10k <- metropolis_hastings(theta0, alpha, beta, 10000)
samples_100k <- metropolis_hastings(theta0, alpha, beta, 100000)

# Histogramas
par(mfrow = c(1, 3))
hist(samples_1k, breaks = 30, probability = TRUE, 
     main = "N = 1000", xlab = "Theta", col = "skyblue")
hist(samples_10k, breaks = 30, probability = TRUE, 
     main = "N = 10000", xlab = "Theta", col = "skyblue")
hist(samples_100k, breaks = 30, probability = TRUE, 
     main = "N = 100000", xlab = "Theta", col = "skyblue")

A medida que aumenta la cantidad de repeticiones, la distribución obtenida se vuelve más estable y representativa de la distribución objetivo. Con pocas iteraciones (N=1000, N=1000), la muestra puede ser ruidosa e imprecisa, con una mayor variabilidad debido a la dependencia inicial de las propuestas y aceptaciones. Sin embargo, al incrementar el número de iteraciones (N= 100000, N=100000), la muestra converge hacia la distribución gamma real, capturando su forma y reduciendo la influencia de las condiciones iniciales.

# Efecto de la condición inicial
samples_theta0_0.1 <- metropolis_hastings(0.1, alpha, beta, 100000)
samples_theta0_10 <- metropolis_hastings(10, alpha, beta, 100000)

# Histogramas
par(mfrow = c(1, 3))
hist(samples_100k, breaks = 30, probability = TRUE, 
     main = "Theta0 = 1", xlab = "Theta", col = "skyblue")
hist(samples_theta0_0.1, breaks = 30, probability = TRUE, 
     main = "Theta0 = 0.1", xlab = "Theta", col = "skyblue")
hist(samples_theta0_10, breaks = 30, probability = TRUE, 
     main = "Theta0 = 10", xlab = "Theta", col = "skyblue")

La condición inicial (\(\theta_0\)) tiene un impacto significativo en las primeras iteraciones del algoritmo, especialmente si \(\theta_0\) está lejos de los valores más probables de la distribución objetivo. Sin embargo, con un número suficientemente grande de repeticiones, el efecto de \(\theta_0\) desaparece gracias a la convergencia de la cadena de Markov. Esto se evidencia en los experimentos realizados: para \(\theta_0 = 0.1\), \(\theta_0 = 1\), y \(\theta_0 = 10\), todas las cadenas convergen a la misma distribución final cuando se ejecutan \(N = 100000\) iteraciones. Por lo tanto, la condición inicial es importante a corto plazo, pero no afecta el resultado final si el algoritmo se ejecuta por suficientes iteraciones.

# Comparación con la distribución gamma real
samples_real <- rgamma(100000, shape = alpha, rate = beta)

# Graficar comparación
par(mfrow = c(1, 1))
hist(samples_100k, breaks = 30, probability = TRUE, 
     main = "Comparación con Gamma real", xlab = "Theta", col = "skyblue")
lines(density(samples_real), col = "red", lwd = 2)
legend("topright", legend = c("MH Samples", "Real Gamma"), 
       col = c("skyblue", "red"), lwd = 2)

Se observa que el histograma obtenido con el algoritmo Metropolis-Hastings para \(N = 100000\) coincide de manera consistente con la densidad de la distribución gamma. Esto valida que la implementación del algoritmo es correcta y que las muestras generadas se distribuyen según la distribución objetivo con un número suficiente de iteraciones.

A work by CC6104